library(tidyverse)
library(here)
library(janitor)
library(raster)
library(sf)
library(tmap)
library(tmaptools)
library(gstat)
gc_dem <- raster(here("data", "gc_dem.tif"))
# Look at this using plot()
plot(gc_dem)
# Let's check the CRS @ = $
gc_dem@crs
## CRS arguments:
## +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# UTM projection is in meters (this kinda sucks)
# Check the extent (bounding box)
gc_dem@extent
## class : Extent
## xmin : 229520
## xmax : 412580
## ymin : 3982160
## ymax : 4100630
# Lets change the projection to lat/long
# Copy the CRS arguments after running gc_dem@crs
# Create a wgs84 with latlong:
wgs_84 <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs"
gc_reproj <- projectRaster(gc_dem, crs = wgs_84, method = "bilinear")
gc_reproj@extent
## class : Extent
## xmin : -114.0417
## xmax : -111.9681
## ymin : 35.94488
## ymax : 37.04945
### Crop raster to a smaller area (bounding box) ******
bounds <- as(extent(-112.4, -112.0, 36.1, 36.3), 'SpatialPolygons')
crs(bounds) <- crs(gc_reproj)
gc_crop <- crop(gc_reproj, bounds)
aggregate exists in baser, and others. Definitely specify that you want the raster:: version.
gc_agg <- raster::aggregate(gc_crop, fact = 10)
plot(gc_agg)
First, convert data to a data frame:
# xy = df makes lat long a separate column in the data frame!
gc_df <- as.data.frame(gc_agg, xy = TRUE)
# ggplot does not, by default, consider projection
# Use coord_quickmap() to quickly augment this
ggplot(data = gc_df, aes(x = x, y = y)) +
geom_raster(aes(fill = gc_dem)) +
coord_quickmap() +
theme_minimal() +
scale_fill_gradientn(colors = c(
"turquoise",
"magenta",
"firebrick",
"white",
"black",
"cyan",
"lightgreen",
"darkred"
))
# First, make a copy of the cropped raster data
gc_hab <- gc_crop
# So, to keep values within a range (1000 - 1500), we gotta make cells outside of that range as NA
# Vertical line = OR
# From gc_hab, anything thats >1500 OR <1000, assign those to NA
gc_hab[gc_hab > 1500 | gc_hab < 1000] <- NA
plot(gc_hab)
tmap_mode("view")
tm_shape(gc_hab) +
tm_raster(legend.show = FALSE,
palette = "plasma")
Read in KS counties shapefile data:
ks_counties <- read_sf(here("data", "ks_counties", "ks_counties_shapefile.shp")) %>%
clean_names()
plot(ks_counties)
# Plot can tell you if things are read in properly
# Now, check the CRS within sf:
st_crs(ks_counties)
## Coordinate Reference System: NA
# Hey! There is none. Let's set one...
st_crs(ks_counties) <- 4326
st_crs(ks_counties)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
# WOOOoooooo.....
plot(ks_counties)
Practice in ggplot()
ggplot(data = ks_counties) +
geom_sf() +
theme_minimal()
Now, let’s read in the rainfall data:
ks_rain <- read_csv(here("data", "ks_rain.csv")) %>%
clean_names()
r doesn’t know that this is spatial data yet. Let’s get it to recognize this as spatial points. Let’s convert it
ks_sf <- st_as_sf(ks_rain, coords = c("lon", "lat"), crs = 4326)
# You can set the crs right in the
plot(ks_sf)
let’s plot this on top of the counties:
ggplot() +
geom_sf(data = ks_counties,
aes(fill = "lightred"),
show.legend = FALSE) +
geom_sf(data = ks_sf,
aes(color = amt,
size = amt),
show.legend = FALSE) +
theme_bw()
careful here…we’re using functions to different packages. gstat doesn’t like sf…so we needa convert a whole bunch of stuff.
ks_sp <- as_Spatial(ks_sf)
# This gives you min endpoints for a grid
bbox(ks_sp)
## min max
## coords.x1 -101.75 -94.63
## coords.x2 37.00 40.00
# Then make a vector with a bunch of points within that grid...
lat <- seq(37, 40, length.out = 200)
lon <- seq(-94.6, -102, length.out = 200)
# Then, turn this into a spatial grid
grid <- expand.grid(lon = lon, lat = lat)
# This is a 200x200 list of pixels. Woo! Now, lets convert it to spatial...
grid_sf <- st_as_sf(grid, coords = c("lon", "lat"), crs = 4326)
# But, kriging functions don't recognize sf class. so...
grid_sp <- as_Spatial(grid_sf)
plot(grid_sf)
ks_vgm <- variogram(amt ~ 1, data = ks_sp)
plot(ks_vgm)
# Guess your nugget! ~ 0.1
# Guess your sill! ~ 1
# Guess your range! ~ 200
ks_vgm_fit <- fit.variogram(ks_vgm, model = vgm(nugget = 0.1, psill = 1, range = 200, model = "Sph"))
# Fit different models! Spherical = Sph, exponential = Exp, gaussian = Gau
plot(ks_vgm, ks_vgm_fit)
ks_vgm_fit
## model psill range
## 1 Nug 0.1021677 0.0000
## 2 Sph 0.9537357 235.1415
ks_krige <- krige(amt ~ 1, ks_sp, grid_sp, model = ks_vgm_fit)
## [using ordinary kriging]
view(ks_krige@data)
# Kriging gives you a predicted amount, plus the variance associated with that prediction
spplot(ks_krige, "var1.pred")
ks_df <- data.frame(ks_krige@data["var1.pred"],
ks_krige@data["var1.var"],
ks_krige@coords) %>%
rename(longitude = coords.x1,
latitude = coords.x2)
# Let's rename these columns...
# Then convert this to an sf object
rain_sf <- st_as_sf(ks_df, coords = c("longitude", "latitude"), crs = 4326)
plot(rain_sf)
ggplot(rain_sf) +
geom_sf(aes(color = var1.pred))
This is no longer kansas! Let’s crop it…
# First read in a geom of kansas...
ks <- read_sf(here("data", "states"), layer = "cb_2017_us_state_20m") %>%
dplyr::select(NAME) %>%
filter(NAME == "Kansas") %>%
st_transform(crs = 4326)
plot(ks)
# Find the intersection of the two...
rain_sf_ks <- st_intersection(rain_sf, ks)
plot(rain_sf_ks)
ggplot(rain_sf_ks) +
geom_sf(aes(color = var1.pred)) +
scale_color_gradientn(colors = c("white","yellow","magenta","purple")) +
theme_minimal()